pro im_process4, path = path, apmed = apmed , filenames = filenames , dark = dark , sky = sky
;+
; NAME:
;     IM_PROCESS4
; PURPOSE:
;     Process NICFPS images via flat-fielding and sky subtraction, then remove non-Gaussian noise.
;
; EXPLANATION:
;
;	Image processing:
;	1. Call up files containing image, median sky, and dark frames
;	2. Create a flat field from the median sky and dark
;	3. Sky subtract and flat field the image
;	4. Remove non-Gaussian noise (two steps):
;		a) Smooth horizontal, high-noise lines
;		b) Plateau sky to zero line-by-line
;	5. Rectify the image by removing distortion in focal plane (NICFPS_UNWARP2)
;		Note that the image distortion applies only to data after
;		the summer shutdown in 2005.
;
;	This program is designed to be used as part of a pipeline, following
;	use of MEDSIGMA_DARKS.pro to create a master dark and MEDSIGMA.pro to create 
;	a median sky, and preceding STACKNCOMBX1 to align the processed images
;	and stack them. 
;
; EXAMPLE:
;     IDL> im_process4
;
; INPUTS:
;	PATH	- gives directory path in which to find pre-processed images
;
; KEYWORDS:
;
;	APMED	- if set, filters for images with _nrm extension (indicates processing
;		by APMED2.pro, where skies have been medianed to mitigate dither residuals)
;
; REQUIRES:
;
;	NICFPS_UNWARP2.pro
;	MMM.pro
;
; MODIFICATION HISTORY:
;	Written by F. Hearty, Univ. of Colorado
;	Added header, comments - K. Willett, Jun 2005
;	Changed MEDIAN to MMM for re-plateauing of rows
;		(better job of finding true sky value in stellar-contaminated rows) - KW, Jun 2007
;	Added PATH, APMED keywords - KW, Jun 2007
;-

; Image files to be processed

if not keyword_set(path) then path = '~/Astronomy/Research/APO/highz/Q4CU04/UT071118/'
if keyword_set(apmed) then begin
	f_filter = '*nrm.fits*' 
	path = '~/Astronomy/Research/APO/night4/'
endif else f_filter = '*.fits'

if ~keyword_set(filenames) then $
filenames = dialog_pickfile(/multiple_files, path = path, filter = f_filter, $
	title = 'FITS files (images) to be processed') $
    else begin
        if stregex(filenames,'.fits') eq -1 then $
            filenames = readcols(filenames,format='(A)')
    endelse
nimages = n_elements(filenames)

; Median sky for images (generated via MEDSIGMA.pro)

;datpath = '~/Astronomy/Research/APO/night4/'
datpath=path
if ~keyword_set(sky) then $
sky = dialog_pickfile(title = 'Median sky FITS file', filter = '*ms*.fits', path = datpath)

; Master dark frame (generated in MEDSIGMA_DARKS.pro)

if ~keyword_set(dark) then $
dark = dialog_pickfile(title = 'Dark FITS file', filter = '*dark*.fits', path = datpath)

; Read in sky and dark, convert to float

sky = double(mrdfits(sky, 0, skyhead, /unsigned))
dark = double(mrdfits(dark, 0, darkhead, /unsigned))

for i = 0, nimages - 1 do begin

   image = double(mrdfits(filenames[i], 0, header, /unsigned))

   ; Subtract sky and flat field

   image_ss = image-sky
   image_ss_ff = image_ss / ((sky-dark)/double(mean(sky-dark)))

   for k=4, 1019 do begin


        ; If line is particularly noisy, then replace it with an average
	; of the lines directly above and below

        if (stdev(image_ss_ff[447:576,k]) gt (5.0 * stdev(image_ss_ff[447:576,k+1]))) $
	then begin
        image_ss_ff[4:1019,k]= 0.5 *((image_ss_ff[4:1019,k-1])+(image_ss_ff[4:1019,k+1]))
        print,'Noisy lines: ',k
	endif

        ; Re-plateau each row

	mmm, image_ss_ff[4:1019,k], row_plateau
        image_ss_ff[4:1019,k]=image_ss_ff[4:1019,k]-row_plateau

  endfor

  ;*****

  ; Rectify the image by removing distortion using NICFPS_UNWARP2

  image_ss_ff=nicfps_unwarp2(image_ss_ff)

  ;*****

  ; Extract filename and strip directory path

   file_nopath = reverse(strsplit(filenames[i], '/', /extract))
   file_begin = strsplit(file_nopath[0], '.fits', /extract, /regex)

  ; Write processed image and original header to file in current directory

   pcr2_filename = file_begin[0] + '_pcr2.fits'
   print, 'WRITING TO ' + pcr2_filename
   mwrfits, image_ss_ff, pcr2_filename, header

endfor

; End program

end
